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t  SUMMARY 

A  general  purpose  ionospheric  ray  uacing  procedure  is  developed.  The  procedure 
is  based  on  a  numerical  solution  to  the  Haselgrove  ray  tracing  equations  and 
includes  calculations  of  ionospherically  induced  Doppler  shift.  The  procedure 
can  handle  a  wide  variety  of  ionospheric  descriptions  and  includes  the  effect  of 
f  the  Earth's  magnetic  field. 


©  COMMONWEALTH  OF  AUSTRALIA  1993 

AUG  93 


APPROVED  FOR  PUBLIC  RELEASE 


POSTAL  ADDRESS:  Director,  Surveillance  Research  Laboratory,  PO  Box  1500,  Salisbury,  South  Australia,  5108 


UNCLASSIFIED 


L 


1 


UNCLASSIFIED 


SRL4)131-TR 


CONTENTS 

Page  No 

1  INTRODUCTION . 1 

2  THEORY . 1 

3  IMK.EMENTATION . 3 

4  THE  PROCEDURE . 5 

5  CONCLUSIONS . 7 

REFERENCES  . 9 

APPENDICES 

A  The  Ray  Tracing  Software . II 


UNCLASSIFIED  iii 


SRL-0131-TR 


UNCLASSIFIED 


Tins  IS  A  BLANK  PAGE 


iv 


UNCLASSfflED 


UNCLASSIFIED 


SRL-0131.tr 


1  INTRODUCTION 

Radiowave  systems  that  rely  on  the  ionosphere  have  a  need  for  reliable  predictions  of  propagation 
performance.  Unfortunately  the  ionosphere  has  an  extremely  complex  motphology  that  is 
.'ontinually  changing  and  so  diese  predictions  are  by  no  means  easy.  The  ionospheric  structure  can 
tc  derived  from  soundings  and  models,  but  the  effect  on  prt^gation  can  only  be  gauged  from 
silutions  to  the  appropriate  electromagnetic  field  equations.  For  HF  radiowaves,  however,  a 
ge  imetric  optics  approximation  is  often  aifequate  and  so  the  propagation  can  be  studied  through  a 
suf  able  ray  tracing  procedure.  Such  a  procedure  needs  to  take  account  of  the  Earth's  magnetic  field 
and  to  accept  a  variety  of  ionospheric  descriptions.  In  addition,  it  wilJ  need  to  calculate  quantities 
such  as  pha^  path,  group  path  and  ionospheric  Doppler  shift  It  is  the  purpose  of  the  current  report 
to  describe  a  computer  procedure  that  can  meet  these  demands. 

Sectirwi  2  of  this  report  presents  the  theoretical  results  upon  which  the  procedure  is  based.  Section  3 
outlines  the  numerical  techniques  that  are  employed  and  Section  4  discusses  the  use  of  the 
procedure.  The  computer  implementation  is  contained  in  Appendix  A. 


2  THEORY 

The  core  of  the  procedure  is  based  on  a  numerical  solution  to  the  Haselgrove  ray  tracing  equations 
<Haselgrove,  1955,  1960  and  1%3).  These  ordinary  differential  equations  (ODE’s)  arc  given  by 

^  =  Ju-KYv,  0) 

dr 


and 


dr 


^£{Ku*MYv) 


d{Yv) 


(2) 


where  X=8.06x10‘-sV/f^,  f  is  the  wave  frequency  {MHz),  N  is  the  electron  density  (electrons 
cm'^).  Y-=fJf  and  f\.  is  the  gyrofirequency  (MHz).  Variable  z  parameterises  the  uajectory,  x  is  a 
position  vector,  u  is  a  vector  that  is  normal  to  the  wavefront  and  /  is  a  unit  vector  in  the  direction 
of  the  magnetic  field.  Qsiantiiies  J,  K,  L  and  M  are  defined  by 

J  ^  2i2(1  -X- Y^p*Y(^  *(v.u)^)  .  K  =  -2X(t/.v5(pl^-1)  (3) 

L  =  .'(1  -X~Y^)p-Y  ,  M  -  2Xp(pV-1) 

where  p  satisfies  the  quadratic 

1  -X-Y^p^*Y{1  *{v.lD^)p-(v.(^^  =  0  (5) 


and  the  different  roots  correspond  to  the  ordinary  and  extraordinary  rays. 
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In  the  current  work,  the  Earth’s  magnetic  field  is  represented  by  a  tiled  dipole 


B  =  -V 


where  Af  is  the  magnetic  moment  (Davies,  1990)  and  is  the  position  of  the  dipole. 
The  phase  path  P  satisfies 

^  ^  ^  oosa  ^ 

<7z  az 


the  group  path  P' 


_3_  =  ti'  cosa _ 

crt  dt 


and  the  ionospheric  doppler  shift  Af  (Bennett,  1%7) 


^  =  -1^  cosa  ± 
ck  c  dt  di 


where  s  is  the  distance  along  the  ray  (ds-  |c^|)  and  a  is  the  angle  between  the  wave  normal  and 
the  ray  direction.  The  refractive  index  n  is  derived  from  the  Appleton-Hartree  equation  (Davies, 
1990) 

>.’■1- - ^  (10) 

2(1  -  X)  -  y*  sityp  ±  i/  r  sin*p  4.  4)^(1-X)*  cos^P 

where  P  is  the  angle  between  the  wave  normal  and  the  magnetic  field  of  Earth,  The  expression  for 

the  group  refractive  index  ji'  is  derived  from  n'=^(fn). 

df 

In  the  case  that  the  magnetic  field  can  be  ignored,  the  trajectories  are  calculated  according  to 


1  dX 

■ST '  "rgr 


where  the  new  x  can  be  interpreted  as  the  group  path. 
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3  IMPLEMENTATION 

The  two  most  important  issues  in  the  implementation  of  the  above  scheme  are  the  choice  of 
ionospheric  representation  and  the  technique  for  solving  the  ODE's.  In  the  present  work,  the  ODE’s 
are  solved  by  a  Runge-Kutta-Fehlberg  (RKF)  scheme  that  can  be  described  as  follows.  Consider  the 
vectOT  system 


(13) 


with  at  time  L.  Form,  in  order,  the  vectors 


(14) 


12 


K  ^  I’ff 


— € 


These  will  yield  fourth  order 


and  fifth  order 


.2, 

'z' 


2!'  =  z  + 


hf{t^*—h  ,  z  +—k  1 
-(*4  4-'J 

(15) 

t,*lh  ,  z  *—k  +  1 

‘  8  -*  32-^  32“* J 

(16) 

1932^.7200.  ^7296^^ 
-*  2197-1  2197“®  2197-13 J 

(17) 

.  ^439^  ^3680^.  845  ^^ 

■k  216"'  513-^  4104"*^ 

(18) 

-®)c.2/r -3544^.  1859^.  11^^ 

27-'  -«  2565"®  4104-*  40"* j 

(19) 

U  .2197^  .  1, 

25^-13  4104-  5- 

(20) 

6656  .  ^  28561  ^  k 

12825*  56435*  55*  55* 

(21) 

Runge  Kutta  approximations  to  the  solution  at  time  f*.,  An  estimate  of  the  error 

incurred  in  the  fourth  order  scheme  is  given  by  the  difference  between  the  fourth  and  fifth  orcter 
approximations.  In  the  RKF  method,  the  steplength  is  adjusted  such  that  this  error  remains  below  a 
preassigned  value.  The  fifth  order  solution  is  used,  however,  for  improved  accuracy.  The  resulting 
scheme  is  extremely  flexible  and  can  adjust  to  the  demands  of  a  highly  varying  ionosphere.  In  the 
present  application,  the  vector  Z  has  nine  components.  These  consist  of  the  position  and  wave 
normal  components  together  with  the  phase  path,  the  group  path  and  the  Doppler  shift.  The 
corresponding  ODE’s  are  given  by  equations  (1),  (2),  (7),  (8)  and  (9). 
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Th/EC  options  for  the  representation  of  ionospheric  electron  density  have  been  incoqxjrstted  into  the 
current  computer  procedure.  The  simplest  of  these  allows  the  electron  density  to  be  defined  through 
the  function  subroutine  ELDEN.  This  is  an  extremely  general  facility  that  can  accept  mostC 
functions.  (Since  numerical  derivatives  arc  used  in  this  option,  it  might  experience  difficulties  if  the 
ionosphere  has  extremely  rapid  variations.) 


The  second  option  provides  for  a  description  in  terms  of  layer  parameters.  These  are  defined  on  a 
regular  geographic  grid  (fixed  increments  in  longitude  and  latitude)  for  two  time  slices.  There  are 
three  layers  (E,  FI  and  F2)  which  are  defined  through  their  heights  (hmE,  hmFl  and  hmF2),  the 
plasma  frequencies  at  these  heights  (foE,  foFl  and  foF2)  and  the  layer  thicknesses  (ymE,  ymFI  and 
ymF2).  The  electron  density  profile  is  composed  of  Chapman  layers,  that  is 


N  =  CM 


2{h-hmE) 

ymE 


1  2{h-hmFl)y^^^ 

2  ymPi  J  " 


1 


•j2(h-hmF2) 


(22) 


where 

Ch(C,x)  =  exp(C(1  -x-exp(-j<)))  (23) 


and  the  coefficients  are  chosen  such  the  electron  density  has  the  correct  values  at  the  layer 
heights  (a  plasma  frequency  of  4  MHz  corresponds  to  electron  density  of  124074*  electrons  per 

cm^).  In  order  to  estimate  the  density  values  between  the  grid  points,  some  son  of  interpolation 
scheme  is  required.  Although  sophisticated  approaches  such  as  cubic  splines  are  desirable,  these  are 
too  inefficient  fex  a  general  purpose  procedure.  Instead,  the  multi-dimensional  interpolation  is  built 
out  of  a  series  of  localised  one  dimensional  cubic  interpolations.  Consider  the  meridian  that  passes 
through  the  point  of  interest  and  calculate  its  intersection  with  the  four  closest  parallels  that  pass 
through  grid  points.  Derive  values  at  the  points  of  intersection  by  use  of  one  dimensional 
interpolations  along  each  of  the  parallels.  The  value  at  the  point  of  interest  is  obtained  by  means  of 
a  one  dimensional  interpolation  between  the  points  of  intersection. 


There  are  several  possibilities  for  the  one  dimensional  interpolation  and  two  options  are  included  in 
the  current  procedure.  The  first  provides  an  approximation  with  C  continuity.  Let  the  grid  points 
he  located  a  unit  distance  t^art  and  the  desired  point  lie  between  the  middle  points  a  distance  u 
from  the  .second  point.  If  the  data  has  values  e, ,  and  at  the  interpolation  poinLs,  the 

value  at  the  desired  point  is  approximately 


2  _  u  _  _  5  O'  u  ^2^2  _  3^,  (24) 

'  rt' 


If  the  data  points  are  sparse,  it  has  been  found  advisable  to  sacrifice  the  C  continuity  and  to  use 
Lagrange  interpolation  instead.  In  this  case,  the  value  at  the  desired  point  is  approximately 


^u{u-1){u -2)  *^{u *^){u -\)iu -2)  u{u -2)  +^{u *^) u{u -^)  (25) 

6  2  2  6 


The  third  option  for  ionospheric  representation  replaces  the  Chapman  layers  with  a  set  of  samples 
above  each  geographical  grid  point.  Sampling  is  uniform  and  so  the  above  geographic  interpolation 
.scheme  can  be  extended  to  include  the  extra  dimension.  This  option  is  by  far  the  most  efficient  and 
is  recommended  when  a  large  number  of  rays  are  to  be  traced. 
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4  THE  PROCEDURE 

The  ray  tracing  procedure  has  been  implemented  as  the  FORTRAN  subroutine  HASEL.  Each  call 
to  this  subroutine  will  trace  rays  that  are  launched  at  a  given  elevation,  in  a  given  direction,  from  a 
given  suining  point  and  with  a  given  wave  frequency.  The  rays  can  be  tra:ed  for  up  to  three  hops 
when  the  effect  of  Earth's  magnetic  field  is  included  and  for  up  to  ten  hops  when  not  A  typical  call 
to  the  subroutine  is  of  the  form 

CALL  HASEL(Long,Lat,Elav,Bear,Freq,Cha,Nos,Out,Type,Tol,M.L) 
where  the  input  parameters  are  defined  as  follows: 

1)  Long  is  a  real *8  variable  that  specifies  the  geographic  longitude  of  the  start  point  (in 
degrees). 

2)  Lat  is  a  real*8  variable  that  specifies  the  geographic  latitude  of  the  start  point  (in  degrees). 

3)  Elav  is  a  real*S  variable  that  specifies  the  initial  elevation  of  the  ray  (in  degrees). 

4)  Bear  is  a  real*8  variable  that  specifies  the  initial  bearing  of  the  ray  (degrees  East  from 
North). 

5)  Freq  is  a  real*8  variable  that  specifies  the  wave  frequency  (in  MHz), 

6)  Tol  is  a  real*8  variable  that  specifies  the  tolerance  for  the  RKF  method  at  each  step  (in 
kilometres). 

7)  Cha  is  a  character  variable  that  has  the  value  'y'  if  the  effect  of  Earth’s  magnetic  field  is  to 
be  included  and  'n'  if  not. 

8/  Nos  is  an  integer  variable  that  specifies  the  number  of  hops  that  are  requned  ithis 
parameter  is  altered  on  exit  from  the  subroutine). 

9)  M  is  an  integer  variable  that  specifics  the  type  of  ionospheric  representation.  The  analytic 
representation  using  subroutine  ELDEN  is  chosen  when  M=7,  the  layer  model  defined  on  a 
geographic  grid  when  Ms6  or  -6  and  the  representation  in  terms  of  height  samples  on  a 
geographic  grid  when  M=10  or  -10.  (Negative  values  of  M  invoke  the  Lagrange 
interpolation  option.) 

10)  L  is  an  integer  variable  that  specifies  the  channel  number  on  which  details  concerning  the 
ray  tfajectories  are  to  be  output  (if  L  <  1  this  information  will  be  omitted).  For  each  step  of 
the  RKF  scheme,  this  option  outputs  the  height  above  the  ground,  the  longitude  ;ind  the 
latitude.  In  addition,  the  angle  between  the  ray  and  Earth's  magnetic  field  is  output. 

The  function  option  (M=7)  defines  the  ionosphere  through  the  real*8  function  subroutine  ELDEN. 
This  subroutine  has  four  real*8  arguments  consisting  of  time  elapsed  (measured  in  seconds  from  the 
start  of  ray  tracing),  longitude  (in  degrees),  lal!*’ide  (in  degrees)  and  distance  from  Earth's  centre  (in 
kilometres).  A  typical  example  of  ELDEN  is  as  follows: 

REAL*8  ELDEN(t,long,lat,r) 
implicit  real*8  (a-l,n-z) 

common  /dat/foE,hmE,ymE,foFI,hmFi,ymFl,foF2,hmF2,ymF2,re,model 

H=l(X).0d0 

hmF2=350.0d0 

foF2=12.0d0 

veN.Ol 

x=(r-6378.  ld(>-hmF2-vel*t)/H 

NF2=12407*foF2**2 

ELDEN=NF2''cxp(  l.()d0-x-exp(-x)) 

return 

end 

This  represents  a  Chapman  layer  ionosphere  with  peak  height  that  rises  at  a  constant  speed. 


UNCLASSIFIED 


5 


S11L-0131.TR 


UNCLASSIFffiD 


The  common  block  is  optional,  but  can  be  used  to  pass  back  infonnation  concerning  ihe  layer 
heights  and  hence  allow  a  meaningful  labelling  of  rays. 

The  grid  options  (M  =  6,-6, 10  or  -10)  require  further  informahon  to  be  passed  from  the  calling 
program  through  a  common  block.  The  following  statement  will  need  to  be  placed  in  the  calling 
program 

COMMON  /grid/MinLaJdinLo,MinR.DLa,DLo.DR,DT,Grid,NR,NLa.NLo 
with  the  parameters  defined  as  follows: 

a)  MinLa  is  a  real*8  variable  that  specifies  the  minimum  latitude  on  the  grid  (ii. 
degrees). 

b)  MinLo  is  a  real*8  variable  that  specifies  the  minimum  longitude  on  the  grid  (in 
degrees). 

c)  MinR  is  a  real*8  variable  that  specifies  the  minimum  distance  (in  kilometres)  of 
the  grid  from  Earth’s  centre  (options  10  or  -10). 

d)  DLa  is  a  real*8  variable  that  specifies  the  increment  in  latitude  on  the  gnd  (in 
degrees). 

e)  DLo  is  a  real*8  variable  that  specifies  the  increment  in  longitude  on  the  grid  (in 
degrees). 

ft  DR  is  a  real*8  variable  that  specifies  the  increment  in  height  on  the  grid  (in 
kilometres). 

g)  DT  is  a  real*8  variable  that  specifies  the  increment  in  time  between  the  time  slices 
(in  seconds). 

h)  GRID  is  a  three  argument  real*4  anay  that  contains  grid  values  of  either  the  layer 
parameters  or  the  electron  density.  The  third  argument  of  the  array  fixes  the  time 
slice  and  can  only  have  the  values  1  or  2.  Geographical  location  is  specified  by  the 
second  argument  As  this  argument  inaeases,  it  runs  eastwards  along  rows  of 
constant  latitude  in  turn.  The  rows  start  at  the  minimum  latitude  and  move 
upwards.  The  meaning  of  the  first  argument  depends  on  the  value  of  M.  For  M=6 
or  -6,  the  argument  values  1  to  9  label  the  layer  parameters  foE,  hmE,  ymE,  foFl. 
hmFl,  ymFl,  foF2,  htnF2,  ymF2  in  order  (f'cqucncies  in  MHz  and  lengths  in 
kilometres).  In  the  case  that  M=10  or  -10,  the  first  argument  labels  the  height 
samples  of  electron  density  (electrons  cm'’)  starting  at  the  lowest  height  and 
moving  upwards.  The  dimensions  of  the  first  ard  second  arguments  should  be  set 
to  the  largest  values  (hat  will  be  required.  (For  ootions  M=6,-6,10  or  -10,  the  array 
GRID  must  be  defined  in  the  calling  program.) 

i)  NR  is  an  integer  variable  that  specifies  the  number  sample  heights. 

j)  NLa  is  an  integer  variable  that  specifies  the  number  of  sample  latitudes. 

k)  NLo  is  an  integer  variable  that  specifics  the  number  of  sample  longitudes. 

The  output  from  HASEL  is  contained  in  the  variables  Out,  Type  and  Nos.  These  are  defined  as 

follows: 

1 )  Nos  is  an  integer  variable  that  specifies  the  total  number  of  rays  U'aced. 

2)  Out  is  a  two  argument  real*8  array  (dimensions  11  and  14  respectively)  that  contains 
infom.ation  about  the  rays.  Each  ray  is  labelled  by  the  second  argument  (values  from  1  to 
Nos).  The  first  argument  labels  the  ray  attributes  as  follows: 

Out(  1.*)  =  Doppler  shift  (Hz) 

Out(2,*)  =  ground  range  (kilometres) 

Out(3,*)  =  group  path  (kilometres) 

Out(4,*)  =  phase  path  (kilometres) 

Ou((5,*)  =  maximum  height  (kilometres) 
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Out(6,*)  =  iniiiai  elevation  (degrees) 

Oui(7,*)  =  final  elevation  (degrees) 

Out(8,*)  =  final  longitude  (degrees) 

Out(9,*)  =  final  latitude  (degrees) 

OutdO,*)  =  bearing  (degrees  East  from  North) 
Outdl,*)  =  wave  frequency  (MHz). 


Type  is  a  character*2l  array  (dimension  14)  that  contains  information  concerning  the  mode 
of  the  ray  that  corresponds  to  the  argument. 


There  is  an  additional  facility  for  situations  where  the  ionosphere  is  periodic  in  longitude.  If  the 
value  of  variable  nos  is  negated,  the  geographic  mesh  (options  m=  10,  -10,  6  or  -6  only)  will  be 
treated  as  one  period  of  a  periodic  structure.  This  facility  is  extremely  useful  for  an  ionosphere  that 
does  not  vary  with  longitude  (set  nlo=l  in  this  case). 


5  CONCLUSIONS 


The  present  report  has  described  a  ray  tracing  procedure  that  can  analyse  the  propagation  of 
radiowaves  through  a  variety  of  ionospheres  and  includes  the  effect  of  the  Earth's  magnetic  field. 
There  are  numerous  uses  for  the  procedure,  mainly  in  the  areas  of  propagation  prediction  and 
ionospheric  research,  in  particular,  two  recent  applications  are  the  simulation  of  backscatier 
ionograms  and  some  investigations  of  the  Doppler  characteristics  of  radiowaves  that  propagate 
through  travelling  ionospheric  disturbances.  A  version  that  estimates  irregularity  induced  Doppler 
spreading  has  also  been  developed  and  this  has  been  used  to  investigate  the  phenomenon  of  Doppler 
spread  clutter. 
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APPENDIX  A  The  Ray  Tracing  Software 


•*•**  HASEL  (version  2.2) 


The  following  routine  performs  lonosherlc  ray  tracing-  It  can  accept 
a  variety  of  Ionospheres  an<3  uses  a  Runge-Kutta-Feblberg  scheme  to 
solve  the  Haselgrove  equations. 


Unless  otherwise  stated,  all  distances  are  In  hllonatres,  angles  In 
degrees,  fre<iuanclea  In  MegaHertz  and  time  In  seconds. 


'**•  input  •••***•*••••**•••**•••••**•*••*•*•••**•••••*•*••••••*•*** 

REAL'S 

along  •  longitude  of  the  start  point  (degrees  B  of  W) 

alat  3  latitude  of  start  point  (degrees,  -ve  for  S  of  Equator) 

elav  s  Initial  elevation  (degrees  from  horizontal) 

bear  ^  initial  bearing  (degrees  east  from  north) 

f  •  wave  frequency 

tol  3  tolerance  (Kms)  at  each  step  of  raytracing  (a  value  of 
l.d-6  Is  sufficient  In  most  cases) 

CHARACTER 

cha  •  'y'  If  magnetic  fields  are  to  be  Included,  'n'  otherwise 
INTEGER 

nos  -  number  of  hops  (changed  on  output  to  total  number  of  rays) 

1  »  channel  number  for  ray  trajectory  Information  (0  If  no 

Information  required) 

m  a  7  for  an  analytic  model  that  Is  described  by  REAL'S 
function  ELDEM 

6  for  a  description  In  terms  : 

of  layer  parameters  on  a  2D  :  Information  passed 
grid  at  two  time  slices  :  via  the  comnon 

10  for  a  description  in  terms  ;  block  GRID 
of  electron  density  values  : 
on  a  3D  grid  at  two  time 
slices  : 


'**•  output  ••••**♦***•**•••••**••*•••*••*•*•***••**•*••****••' 

REAL'S 

out(i,*)  »  doppler  shift  (Hz)  : 

out (2,*)  •  ground  range  : 

out  ( 3 ,  • )  =  group  path  •. 

out (4,*)  a  phase  path,  : 

out (5,*)  a  maximum  height  : 

out(6,*)  a  Initial  elevation  :  '  refers  to  the  ray  number 

out{7,')  a  final  elevation  : 

out (8,*)  a  final  longitude  : 

out (9,*)  a  final  latitude  : 

out (10,*)  a  bearing  : 

out (11,*)  a  frequency 

CHARACTER*21  : 

typo  ( • )  a  mode 

INTEGER 

nos  a  number  of  rays  traced 

Mote: 

The  program  outputs  the  ray  trajectory  and  values  of  the 
angle  between  the  ray  and  Earth's  magnetic  field.  Bach 
line  of  output  will  consist  of  a  height,  a  longitude,  a 
latitude  and  an  angle. 


**  subroutines  required 

LAYERHT,  RKF,  FUNC,  LAYERX,  TBRP,  DTERP,  ELDEM,  ELECTRONR,  ELOEMX 
ELECTRON,  MAQfBTIC 


*"*  common  bloc)c  (  to  be  set  by  user  under  certain  options) 
The  parameters  of  the  common  block  GRID  are.  In  order, 
REAL'S 

blhla  a  lowest  value  of  latitude 
blhlo  a  lowest  value  of  longitude 

rmg  •  height  of  Ionosphere  above  the  Barth  centre 
dla  a  latitude  Increment 

dlo  «  longitude  increment 
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drg  a  haight  l&crament  ({or  s  •  10  only) 
dtz  •  tlaa  batwsan  temporal  ellcee 
INTEGER 

nr  «  niuBber  o£  halgbt  aaaples  for  m  >  10  or  9  for  m  •  6 
nla  a  number  of  latitude  samples 
nlo  =  numhar  of  longitude  samples 
REAL* 4 

yp(ll,12,13}  a  electron  density  samples  if  m  a  lo 
or  layer  parameter  samples  If  m  a  6 

for  ma6  yp  contains  the  ll*th  layer  parameter  In  the  list 

foE,bmE,ymE,foFl,fcmFl,ymFl. f oP2 , hmF2 , ymF3  (frequencies 
In  MHz  eutd  heights  In  Km) 

for  maio  yp  contains  electron  density  (electrons  per  cubic  cm) 
samples  at  height  rmg+ (11-1) “dro 

12  a  (l-l)*nlO'fj  for  the  samples  at  latitude  blhla'f (l-l)*dla 
and  longitude  blblo+ ( j-l) *dlo 

13  labels  the  time  slice  (1  or  2) 

Notes : 

1}  Maximum  values  of  11  and  12  have  been  set  at  61  and  961.  If 
these  limits  are  to  be  changed,  remember  to  alter  for  all 
occurences  of  yp. 

2)  A  maximum  value  of  lOXm  Is  advised  for  drg  ,  S  degrees  for  dla 
and  5  degrees  for  dlo. 


switches  •*•*•••••***••*••••••••••*••*•**••*••*••*«•*•••••**•( 

1)  If  tol  -  -tolerence  the  program  will  execute  faster,  but  the 
doppler  will  be  baaed  on  a  stationary  Ionosphere. 

2)  If  1<1  no  ray  path  data  will  lae  output. 

3)  If  m<0  Iiagrange  Interpolation  Is  used. 

4)  If  noa<0  the  Ionosphere  will  be  given  a  periodic  extension  In 
longitude. 


common  bloOcs  that  are  internal  *••••**•*••*•****•••*••****•**•• 

DAT  holds  (In  order)  current  values  (real*8)  of  layer  parameters 
foE,  hstS,  ymE,  foFl,  hmPl,  ymPl,  foF2,  bmF2  and  ymP2.  The  final  two 
parameters  are  the  Earths  radius  (real*e  In  Kms)  and  the  current 
model  (an  Integer) . 

APPROX  consists  of  a  single  character  variable  that  has  value  'L' 
for  Lagrange  Interpolation  and  'S’  for  smooth  jointed  Interpolation. 
MAG  has  eight  real*8  variables  followed  by  an  Integer  variable.  The 
real  variables  are  the  current  gyrofrequency  (MHz),  the  wave 
frequency  (MHz),  the  Cartesian  coordinates  (Earth  centred)  of  the 
magnetic  dipole  and  the  dipole  moments.  The  integer  variable  has 
the  value  0  when  the  magnetic  Is  to  be  Ignored,  a  value  1  when 
ordinary  rays  are  to  be  traced  euad  value  -1  when  extraordinary  rays 
are  to  be  traced. 


a ubrout Ine  KASEL ( along, alat , elav, bear, f , cba , nos , out , typo, tol , m, 1 ) 
Implicit  real*8  (a-h,o-z) 
roal*4  yp(61,961,2) 

character  type*2l, tYpo*21,chl,cH*21,cl2*21,cha,app 
dimension  xdl)  ,ytmp(ll)  ,xx(ll)  ,xtmp(ll ) ,  out  (11, 14 ) ,  typo (14) 
6,v(3),ryx(3,3) 

common  /dat /£oE,hmE, ymE, foFl,hmFl,ymFl,foF2,hmF2,ymF2, re, model 
common  /grld/blhla,blhlo,rmg,dla,dlo,drg,dtz,yp, nr, nla, nlo 
common  /approx/app 

common  /mag/fhs, £r,xm, ym, zm, px,py, pz, ir 
re»6378.13S 
scale«4S./atan(l. ) 
nv9 

noz«noa 

fr-f 

l£(m.lt.O)  then 
app» ' S ' 
else 
app»'L' 
endlf 

modal-aba (m) 
l£(miodel.eq.7)  dtz-1. 
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l£(tol.lt.0.0)  then 
dtwadtz 
dtZaO.O 
endlf 

tolx=abs(tol) 

*****  dipole  moment  «  (pxt,pyt,pzt)  and  dipole  location  =  (xin,ym,  zm)  *** 
dlat=77.75 
dlona-295.75 

pxtecoa (dlat/acale} *cos (dlong/scale) 

pytscos (dlat/acale) *sln(dlong/Bcale) 

pzt^aln (dlat/acale) 

xm— 434.66 

ymsl99 .47 

zm«BO .25 

fha»2.8*.31*re**3 
pxa£ba*pxt/£r 
py»  f ha ‘pyt / £r 
pz=£ha*pzt/£r 
1£ (cha.ne. 'y' )  £h8>0.0 


tzhla-blhla+real (nla-1) 'dla 
trhlo^blhlo+real (nlo-1) “dlo 
nhopaab8(noa) 

If (cha.eq. 'n' )  then 

lra=0 

lrb=0 

elae 

lra=l 

lrbs2 

l£(ab8 (noa) .eq. 2)  nhop>3 
If (aba (noa) .qe. 3 )  nhop*7 
endlf 
noaaO 


do  909  lrr>lra, Irb 

*****  ray  type  loop  *•*••••*••••*•*****•***•*•*••••*•••*•*•*•••*•••• 

If (Irr. eq.O)  then 
ohl«'N' 
lr»0 

elae  l£(lrr.eq.l)  then 

ohl='0' 

lr=l 

else  l£(lrr.eq.2)  then 

chl='X' 

lr«-l 

endlf 

drseln(elav/scale) 

dt--cos (elav/acale) *coa (bear/scale) 

dp>coB(elav/acale)*8ln(bear/scale) 

phlsatand.  )*along/45. 

tbeta>2 .*atan(l . ) -atan(l. )*alat/4S. 

x(l)a:ra*8ln(theta)*co8(pbl) 

x(2)zre*aln( theta )*sin(phl) 

x(3)are*cos (theta) 

xl=x(l) 

x2-z(2) 

x3=x(3> 

x(4)«dr*aln( theta ) *  coa ( phi ) tdt  *co8 ( theta ) *co8 ( phi ) -dp*  a in ( phi ) 

x(5)>dr*8ln(tbeta)*8in(phl)  'fdt*cos(tbata)  *aln(phl)  'fdp*coa  (phi) 

x(6) «dr*coB (theta) -dt*aln(theta) 

t-0.0 

lta>2 

l8tep>l 

type- ' 

X(7)-0.0 

x(6)-0.0 

X(9)-0.0 


do  979  lhop«l,ahop 

hop  loop  *•••*••***•**••••**•*** 

zmax-0 . 0 

8c-8qrt(x(4)**2+x(5)**2tx(6)**2) 

x(4)»x(4) /ac 

x(5)«x(5)/bc 
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x(S)«X(6)/8C 

If  (cba.«q. 'y  .amd.lhop.gt.l)  than 
If (lh.op.e<i.2)  than 
do  k>l,n 
xtaip(k)  •xdc) 
enddo 
Ir.l 
chl-‘0' 
else 

If ( Ihop . aq , 3 . or . lhop.eq.4 .or . Ihop.aq. 6)  lr>-l 
if ( ihop.eq. 5 .or . lhop.aq.7 )  Iral 
If ( Ihop . aq . 3 )  than 
do  kal,n 
ytinp(k)»x(k) 
x(  k)axtitip(k) 
enddo 
cllatype 
Ital-ltb-fl 
endlf 

If (Ihop.eq. 4)  then 
do  kal,n 
xtmp ( k 1 ax ( k ) 
enddo 
cllatype 
ItaJaltb*! 
endlf 

Ifdr.aq.l)  chla'O' 
if(ir.oq.-l)  ohl-'X' 

If (ihop.eq. 4 .or. Ihop.eq. 5)  then 
do  kal,n 
X(k)aytBip(k) 
enddo 
typoacll 
Itaaltal 
endlf 

If (ihop.eq. 6 .or. Ihop.eq. 7)  then 
do  kal,n 
x(k)axtsip(k) 
enddo 
typaacl2 
ltaalta2 
endlf 
endlf 
endlf 
zad=0 .0 
zedaO  .0 
Itbalta+l 

type( Ita: itb+1 ) a ' g ' //chi// '  ' 
ltx=2 

If (chi .eq. 'N' )  then 

type(ltb:ltb)a'  ' 
ltb=itb-l 
ltxaltx-1 
endlf 

•••*•  start  of  calculations  for  an  Individual  hop  *••••••••••••••’ 

999  Istepalsteptl 
do  lal,6 
XX(1).X(1) 
enddo 
pi-x(7) 

dOpXaX(8) 

pp-x(9) 

zbdazad 

zadazed 

xlongarlong 

xlatarlat 

ODB  solver  aaa»aa»aaa»aa»a»aaaaaaaaaaaaa»aaaaa«aaaaaaaaaaaa 

hmaxaSO . 
bmlna .01 
hahmax 

call  rkt(t,x,n,h,tolx,hnlp, hniax) 

*  xd)  to  x(3)  a  Cartesian  coordinates  of  currant  ray  position 

*  x>7i  to  x(9)  a  group  path,  dopplsr  shift  and  phase  path  at  this 

*  location 
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red-aqrt  (x(l)  •*2+x(2)  ••2+jc(3.  **2) 

zedlszed 

zadsred-re 

rlat=90 .-acoa (x{3 ) /red) ‘scale 
if  (aba(x(l))  .lt.l.Od-10)  the.r 
rlon(r»atan2  (x(2 ) ,  1.  Od-lO) ‘scale 
else 

rlon3aatan2(x(2),x(l) )‘acale 
endlf 

if  (rlong.  It  .0 .0 )  rlong>rlong4  3(;0. 0 
Ifd.gt.O)  then 

call  inagiietlc(x,v,ryx,ba,dlF.dcl,n} 
vu=x(4)*v(l)+x(5)‘v(2>+x(6)*w(3> 
xta8qrt(abs(x(4)“2+x(S)“2  +  i£(6)“2)) 
angleoacos( . 9999999‘vu/ (ha‘xt ) ) ‘scale 

If (zed. gt. 0.0)  write ( 1, fmts '  4fl5.7)')  zed, rlong, rlat, angle 
endlf 

If ( (rlat.gt.trhla.or.rlat.lt.blnla) .and.model.ne.7)  go  to  909 
If ( (rlong . gt . trhlo . or . rlong . It . nlhlo) 

& .emd.model .ne .7 .and.noz .gt .0)  gc  to  909 


““*  At  the  current  stage  of  solution 
zad=  altitude  of  the  ray 
“*“  rlong  »  longitude  of  ray 
rlat  =  latitude  of  ray 


label  for  highest  reflection  layer  of  current  hop 
If  (zad.gt.max(zbd,zed,zznax) )  tber 
If (model. eg. 10}  then 

call  layerpana(xlong,xlat,hl,lL2,h3, fl,f2, f3) 
bmE^hl 
hmFl=h2 
bmF2:>h3 
endlf 

If (zad. le.hmE)  then 
Itbalta+l 

type(lta:ltb+l)»'E'//chl//'  ' 
ltX32 

else  If (zad.le.hmFl)  then 
ltb«lta'f2 

type (itailtb)='Fl‘/ /chi 
ltx=3 

else  If  (zad.  le.haiF2)  then 
ltbslta+2 

typo ( ita ; Itb) » ' P2 ' / /chi 
ltx»3 
else 

Itbalta+l 

type(ltajltb*l)»'S'//chl//'  ' 
ltx»2 
endlf 

If (chi .eg. 'N' )  then 
type (Itb: Itb) »'  ' 
ltb.ltb-1 
ltx>itx-l 
endlf 
endlf 


zmax«max (zad, zmax ) 

If (zed. go. 600. )  go  to  909 
If ( zed.gt .0. )  go  to  999 
If  ( Ihop.eq.  1  .or  .cba.ne.  'y' )  Itaxlta-flts 
Interpolation  for  values  at  ground  level 
popozedl/ (zedl-zed) 
do  1 > 1 , 6 

x(i)-xx(l)+pop‘(x(l)-xx(l)) 
enddo 
xl-0.0 
xxl>o .0 
cosa-O.O 
do  1-1,3 
xl-xl+x(i)“a 
xxl-xxl+x(l+3)“2 


UNCLASSIFIED 


IS 


SRL-0131-TR 


UNCLASSIFIED 


co8a«cosa+x(l)*x{l+3) 

enddo 

elav£-abs(90.-45.*acoa(coaa/sqrt(xl*xxl) ) /atand. ) ) 

sd=8qrt ( (x(l) -xl)»*2+ (x (2 ) -x2 ) ••2+ (x(3) -x3 ) *»2 ) 

sd-2 . *re*asln( .5*8d/re) 

pl=pl+pop*(x(7) -pi) 

pp=pp+pop* (x ( 9 ) -pp) 

dopx=dopx+pop* ( X ( 8 ) -dopx ) 

x(7)=pl 

X ( 8 ) =dopx 

x<9)»pp 

dopxsl . 0d6*dopx 

xlongsxlong+pop* ( rlong-xlong) 

xlat>xlat  +pop* ( r lat -xlat ) 

no8sno8+l 

out ( 1 , no8 ) >dopx 

out (2, nos) sad 

out ( 3 , nos ) spl 

out (4.ao8)spp 

out (5, nos ) szmax 

out(6,no8)selav 

out ( 7 , nos ) sslavC 

out  <  8 , nos ) sxlong 

out ( 9 , nos ) sxlat 

out  <10, nos) shear 

out(ll,nos)sfr 

typo (nos) -type 

ss-x(l)**2+x(2)**2+x(3)**2 

xr=x(4)»x(l)-fx(5)*x(2)+x(6)*x(3) 

x(4)»x(4)-2.*xr*x(l)/88 

x(5)»x(5)-2.*xr*x(2)/a6 

x(6)=x(6)-2.»xr*x(3)/88 

If(l.gt.O)  then 

zod-0 . 0 

write(l,fnit=' (4fl5.7) ' )  zod,xlong, xlat, angle 

endlf 

979  continue 

909  continue 

If (tol. It .0 . 0)  dtz-dtw 

return 

end 


LAYERFAHH 


layer  parameters  at  a  given  location 


’•••  Input  ••••••••••••••••••••••••••••• 

roal*8 

elong  -  longitude  (deg)  of  sample  point 
elat  -  latitude  (deg)  of  sample  point 
'*•*  output 
real's 

hi  -  height  of  E  layer 

h2  -  height  of  FI  layer 

h3  -  height  of  F2  layer 

fl  -  critical  frequency  of  E  layer 

f2  -  critical  frequency  of  PI  layer 

f3  -  critical  frequency  of  P2  layer 


subroutines  required 
BLDEKX 


common  blocb  **♦*•* 
GRID  (defined  In  HA8EL) 


subrout Ine  layerparm ( elong, elat , hi , b2 , h3 , f 1 , f 2 , f 3 ) 

Implicit  real's  (a-h,o-z) 
real's  yp(61,981,2) 

conmon  /grld/blhla,blhlo, Erag,dla,dlo,drg,dtz,yp,nr,nla,nlo 
re-6378.135 
rO-zmg 

eO-eldenx (elong, elat, rO,dlO,dthet,dpbl, 1) 
rl-rO-fdrg 

el-oldenx<elong,elat,rl,dll,dthet,dphl  1) 
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hO>iisg-re 
hl-0.0 
h2>0.0 
h3-0.0 
ea«0 . 0 
ebsO . 0 
ec>0.0 
do  1=3, nr 
r2=rina+dr3*  ( i-1 ) 

e2=eldenx(elong,elat, r2.dl2,dtbat.dpbl, 1) 

d21=(dl2-dl0) /(r2-r0) 

bb-0.0 

ee=0.0 

bz=0. 0 

if (ail*dl0.1a.l.d-33,and.ab«(dll-dl0) .gt.l.d-33)  then 
checlc  for  waxlmiim  *••••**•***•*•*****•*•****•*••*•••*•••••••• 

hb= (dll»rO-dlO*rl ) / (dll-dlO ) -re 
ea» (dll*eO-dlO*al) / (dll-dIO ) 

If (hh.gt. 100.0)  then 
if {hh.lt. 140.)  then 
hl=hh 
ea>ee 
else 

if (hh.lt. 230. )  then 
h2=hh 
eb=ee 
else 
h3=hh 
ec=ee 
endif 
endlf 
endif 
else 

cbec)c  for  point  of  inflection  •**•**•*********•*•**••*••••*•* 

if(i.gt.3)  then 

if (rl-re.lt.230. )  then 

if (d21*d20 . la . 1 .d-33 .and.abs (d21-d20 ) .gt.l.d-33)  than 
hz= (d21*r0-d20*rl ) / (d21-a20 ) -re 
ez= (d21»e0-d20*el) / (d21-d20) 
if (hz.gt. 100.0. and. hz.lt. 230.)  then 
if (hz.ge.l40. )  then 

lf(h2.1t.l.)  than 
h2=hz 
eb=ez 

endif 

else 

if (hi. It. 1.)  then 
hl>hz 
ea«ez 

endif 

endif 

endlf 

endif 

endlf 

endlf 

endlf 

dl0=dll 

dll-dl2 

d20-d21 

e0>el 

el=e2 

r0«rl 

rl«r2 

enddo 

if(hl.lt.hO)  hl>h0 
If(h2.1t.b0}  b2-bl 
H(h3.1t.h0)  h3»h2 
lf(ec.lt.l.)  then 
ac«e2 
h3=r2-re 
endlf 

fl=sqrt(ab8(ea*80.6d-6) ) 
f2-sqrt(aba(eb*80.6d-6) ) 
f 3=sqrt (abs(ec*80 .6d-6) ) 
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return 

end 

ELECTRON 


electron  density  and  Its  gradient 


*****  input  •*•*•***•**••••**•*•••••••••••*•••****•**••***•••* 

•  real's 

•  x(')  3  Cartesian  coordinates  (Earth  centred)  of  sample  point 
*****  output  **•***********•**••***•••**•••*••••*•*•*•••••••** 

•  real'8 

•  en  =  electron  density  (electrons  per  cubic  cm) 

•  dnx(*)  =  Cartesian  derivatives  of  electron  density 


*****  sxJsroutlnes  re<iuired  *•••••♦••• 
*  ELOENX,  BLDEK,  LAYERX  and  ELECTRONR 


*****  c<xnmon  blocks  ************* 

•  GRID  and  DAT  (defined  in  HASEL] 


subroutine  electron (x,en,dnx) 

Implicit  real's  (a-h,o-z) 
dimension  x(3).dnx(4) 
real'4  yp(61,961,2) 

conmon  /grld/blhla,blhlo< rmg,dla.dlo,drg,dtz,yp, nr, nla.nlo 
comaon  /dat/foE, hmE.ymE, foFl.hmPl. yniFl, foF2, hmF2 ,ymF2, re, model 
r=8qrt(x(l)"2+x(2)"2*x(3)"2) 
scales4S.  /atand . ) 
theta=aco8(x(3) /r) 
phl*atan2 (x(2 ) , x(l) ) 
elat*90 . -scale'theta 
elong*scale'phl 

If ( elong.lt .0.0)  elong*elong*360 . 

*****  Parameter  da  should  be  less  than  20%  of  the  smallest  geographic  " 
*****  scale  (In  degrees).  *• 

da*. 002 

If (model. eq. 10)  then 

en*eldenx(elong,elat, r,dlffr,dthet,dphl,l) 
dthet*scale'dthet 
dphl*scale'dphl 
if (dtz.gt.l.d-20)  then 
enp*eldenx( along, elat, r,dl,d2,d3, >2) 
enm=en 
endlf 
else 

If (model .eq. 7)  then 

*****  Parameter  dr  should  be  less  than  20%  of  the  smallest  height  '**•" 
*****  scale  (in  (cilometers) .  ****** 

dr«  .2 
t=0.0 

snpp*elden(t, elong*da,elat,  r) 
enpm*elden( t , elong-da, elat , r ) 
entp*elden(t, along, elat tda, r) 
entm»elden(t, along, elat -da, r) 

If (dtz .gt . 1 .d-20)  then 
enp*elden (t* . SdO'dtz, along,  elat ,  r) 
enmaelden(t- . SdO'dtz, along, elat, r ) 
endlf 

dlffr* .5*(elden(t, along, elat, r*dr) -eldan(t,elono,alat, r-dr) ) /dr 
else 

call  layerx(elong*da,elat, 1) 
call  electronr(r,enpp,dpp) 
call  layerx(elong-da,elat,l) 
call  electronr(r,enpm,dpai} 
call  layarx(elong,elat*da,l) 
call  elactronr(r,antp,dtp) 
call  layerx(elong,alat-da, 1 ) 
call  alactronr (r, sntm,dtm) 
if (dtz .gt . 1 .d-30)  than 
call  layerx(elong,elat, -2) 
call  alactronr (r,enp,dq) 
anm> . 25' (anpp+anpm*antp*antm) 
endlf 
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dlf£r».25*(dpp-i  Ipm-fdtp+dtm) 
endif 

ens  .25*  (eapp-i-enpm-fentp+entm) 
dphl» .5*8cale* (enpp-enpm) /da 
dthets- .5* scale* (entp-antm) /da 
endif 

if (dtz.lt. l.d-20)  then 
dtiiBe=dphl/  (240  .*  scale) 
else 

dtlmea (enp-enm) /dtz 
endif 

dna(l)= (dlf fr*x(l) +dthet*cos (theta) *co8 (phi ) 
&-dphl*8in(phl) /aln(thata) )/r 
dnx(2)* (dlf f r*x(2) +dthet*coa(theta) *8in(phi ) 
4+dphl*coa (phi) /aln(theta) ) /r 
dnx (3 ) » (dlf £r*x ( 3 ) -dthet*8in( theta) ) /r 
dnx(4)=dtline 
return 
end 

LAYERX  ************** *********************** 


layer  parameters  from  grid  values 


*****  input  **••*•*•**•••***••*••••••**** 

*  real*S 

*  along  =  longitude  (deg)  of  sample  point 

*  elat  s  latitude  (deg)  of  sample  point 

*  Integer 

*  luu  s  time  slice  label  (1  or  2) 


•••••  subroutine  required 
*  TEHP 


*••••  common  bloc)cs  ************* 

*  GRID  and  DAT  (defined  In  HASEL) 


*****  note  •****•*****•***********************' 

*  The  output  la  placed  In  the  ccanmon  bloc]c  DAT 


sxibroutlne  layerx(elong, elat,  luu) 

Implicit  raal*8  (a-h,o-z) 
real*4  yp(61,961,2) 
dimension  yz(4),yy(9) 

common  /grld/blhla, blblo, rmg, dla, dlo, drg, dtz , yp, nr , nla,nlo 
comaon  /dat/yy, re, model 
llp*ab8(luu) 
dlongxelong-blhlo 
If (dlong.gt .360 . )  dlongsdlong-360 . 
la*mln(max(2, Int ( (elat-blbla) /dla) *1) ,nla-2 ) 
loslnt (dlong/dlo) +1 
lo»lo-int( (lo-l)/nlo)*nlo 
If(lo.lt.l)  lolo-fnlo 
deLxs(elat-blhla) /dla-real (la-1) 
dox*dlong/dlo-real ( lo-l ) 
nOs ( la-2 ) *nlo 
nl«?iO+nlo 
n2»nl+nlo 
n3*n2*nlo 
do  1-1,9 
do  )c-0,3 
lt)c-lo-l*)c 

If  OcA.  It .  1)  kk-AA*nlo 
If (kk.gt.nlo)  kk-kk-nlo 
xO-yp ( 1 , n0*kk, 1 Ip) 
xl-yp(l,nl'fkk,  lip) 
x2-yp(l,n3+kk, lip) 
x3-yp(l,n3*kk, lip) 
yz ( k+ 1 ) -terp ( dax, xO , xl , x2 , x3 ) 
enddo 

yy(i)-terp(dox,yz(l) ,yz (2) ,y*(3 ) ,yz (4 ) ) 
enddo 
return 
end 
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ELOENX 


*****  electron  density  and  gradients  from  grid  values 


Input  ***************************** 

*  real*8 

*  elong  =  longitude  (deg)  of  sanvle  point 

*  elat  «  latitude  (deg)  of  sample  point 

*  r  =  distance  from  Earth  centre  (Kms) 

*  Integer 

*  luu  s  time  slice  label  (1  or  2) 


*****  output 

*  real's 

*  eldenx  •  electron  density  (electrons  per  cubic  cm) 

*  diffr  >  derivative  In  vertical  direction 

*  dthet  3  derivative  In  counter  latitudinal  direction 

*  dphl  =  derivative  In  longitudinal  direction 


•••**  subroutines  required 
•  TERP,  DTERP 


*•*'*  common  bloc)c  ****** 
•  GRID  (defined  In  HASEL) 


real's  function  eldenx(elong, elat, r, diffr. dthet, dphl, luu) 
implicit  real's  (a-h, o-z) 
real'4  yp(61,961,2) 

dimension  yz (4) ,yzt (4 ) , e(4 ) ,dt (4) ,dp(4) 
conmon  /grld/blhla,blhlo, rmg,<31e,dlo,drg,dtz,yp,nr,nla,nlo 
llp=abs (luu) 

1  f  ( r .  1 1 .  rmg .  or .  r .  g  t .  rmg  '(nr-D'drg)  then 
eldenx«0. 0 
diffr. 0.0 
dthet. 0.0 
dphl.O. 0 
else 

dlong.elong-blhlo 

If (dlong.gt.360. )  dlong.dlong-360 . 

la.min(max(2, lnt( (elat-blhla) /dla) '1) .nla-2 ) 

lo.lat(dlong/dlo) +1 

lo.lo-lnt ( ( lo-l ) /nlo) 'nlo 

If(lo.lt.l)  lo»lo+nlo 

lr.mln(max(l, int( (r-rmg) /drg) +1) ,nr-2) 
dax. (elat-blhla) /dla-real (la-1) 
dox.dlong/dlo-real ( lo-l ) 
drx. (r-rmg) /drg-real(lr-l) 
n0= (la-2) *nlo 
nl.nO+nlo 
n2.nl 'nlo 
n3.n2'nlo 
lf(lr.lt.2)  then 
nstart.O 
else 

nstart.-l 

endlf 

do  l.nstart,2 
do  )c.0,3 
)ck.lo-l']( 

IfOOt.lt.l)  kk-kJitnlo 
if (kk.gt.nlo)  kk«)tk-nlo 
xO.ypdr'l.nO'kk,  lip) 
xl«yp(lr+l,nl'kk, lip) 
x2.yp(lr'l,n3'kk, lip) 
x3.yp(lr'i,n3'kk, lip) 

1 f ( luu .gt.O)  yzt(k'l) .dterp (dax, xO , xl , x2 , x3 ) 
yz(k'l) •terp(dax, x0,xl,x2,x3) 
enddo 

yy-terp(dox,yz(l) ,yz(2) ,yz(3) ,yz(4) ) 

If (luu. gt.O)  then 

yyt .terp ( dox, yzt(l),yzt(2),yzt(3), yzt  <  4 ) ) 
yyp«dterp(dox,yz(l) ,yz(2) ,yz(3) ,yt(4) ) 

else 

yyt.O. 
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yyp-0. 
endlf 
e(l+2)«yy 
dt ( i+2) =yyt 
dp(l+2)»yyp 
enddo 

if (lr.lt. 2)  then 
e(l)=0.0 
dt(l)=0.0 
dp(l)=0.0 
endlt 

eldenx>terp(drx,e(l) ,e(2),e(3),e(4)) 

If (eldenx.lt. 0.0)  then 
eldenx>0.0 
dlf£r«0.0 
dthetaO.  0 
dphl>0.0 
else 

Ifduu.gt.O)  then 

dl£fr=dterp(drx, e(l),e(2),e(3),e(4)) 
dthet--terp(drx,dt (1) ,dt(2),dt(3),dt(4) ) 
dphlsterp (drx, dp ( 1 ) , dp (2 ) , dp (3 ) , dp (4 ) ) 
else 

dlffr=0.0 
dthetsO . 0 
dphl^O .0 
endlf 
endlf 
endlf 

dthet>dthet /dla 

dphlsdphl/dlo 

dl££r»dlffr/dro 

return 

end 

**•*  TERP  •••**♦****•****•••••****••****•••**•••< 


cubic  Interpolation 


'•••  Input  ***••••••••••••••••••••*•••*••*••*•******•*•***••*♦* 

real*8 

pO,  pi,  p2  and  p3  =  samples  (In  order)  at  points  spaced  a  unit 
distance  apart. 

d2  >  distance  of  teat  point  from  second  sample  p-olnt. 


•••  output  •••••••••*•*•*• 

real's 

terp  «  value  at  test  point 


'•••  common  block 

APPROX  consists  of  a  single  character  variable  that  has  value  'L'  * 

for  Lagrange  Interpolation  and  'S’  for  smooth  jointed  Interpolation.  * 


real's  function  terp(d2.p0,pl,p2,p3) 

Implicit  real's  (a-h,o-z) 
character  app 
coimnon  /approx/app 
If (app.eq. 'L' )  then 
dl-d2'fl. 
d3-d2-l. 
d4-d2-2. 

terp- (P3'dl'd2'd3-p0'd2'd3'd4 ) /6 .  +  (pl'dl'd3'd4-p2'dl'd2'd4 ) /2 . 

else 

terp-pl+ .  5'  (p2-p0 )  'd2  +  (pO-2 . 5'pl+2  .  'p2-  .  5*p3 )  'd2"2 
k*  { .  5'p3-l .  S'p2  +  1 . 5'pl-  .  5'pO  )'d2"3 
endlf 
return 
end 
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•••*•  DTBRP 


•  darlvatlv*  of  cubic  lntari>olacloii 


*****  input  *••*•••***•****»*•****•**•*»*•* ********************* 

*  real* 8 

*  pO,  pi,  p2  and  p3  «  saaplaa  (in  order)  at  points  spaced  a  unit 

*  distance  apart. 

*  d2  *  distance  of  test  point  from  second  sample  point. 


*****  output  ********************* 

•  real's 

•  dterp  «  derivative  at  test  point 


*****  common  block  ***************************************************** 

*  APPROX  consists  of  single  character  variable  that  has  value  'L'  * 

•  for  Lagrange  interpolation  and  'S'  for  smooth  jointed  interpolation.  * 

real's  function  dterp {d2 , pO .pi, p2,p3 ) 

Implicit  real's  (a-h, o-z) 
character  app 
common  /approx/app 
If {app.eq. 'L' )  then 
dl.d2+l . 
d3*d2-l. 
d4=d2-2 . 

dterp- {p3' {d2*d3+dl*d3 +dl*d2 ) -pO' (d3*d4 +d2*d4 +d2*d3 ) ) /6 . 

(pi* (d3*d4+dl*d4+dl*d3 ) -p2'(d2*d4*dl'd4+dl*d2) ) /2 . 
else 

dterp- . 5* (p2-p0 ) + (2 . *p0-5 . 'pi +4 . *p2-p3 ) *d2 
t+{l  .5*p3-4 . 5*p2f4 .5'pl-1.5*p0)*d2"2 
endlf 
return 
end 

*****  det 


30  determinant 


*****  input  ***********•**•"*' 

*  real's 

*  al(*)  -  1st  column  of  matrix 

*  a2(*)  -  2nd  column  of  matrix 

*  a3(*)  »  3rd  column  of  matrix 


*****  output  •*"" 

*  real's 

*  det  -  determinant 


real's  function  det(al,a2,a3) 

Implicit  real's  (a-h,o-z) 
dimension  al (3 ) .a2 ( 3 ) . a3 ( 3 ) 

det-al(l)'(a2(2)*a3(3)-a3(2)*a2(3))-a2(l)*(al(2)*a3(3) 
i-al(3)*a3{2) )+a3(l)*(al(2)*a2{3)-al{3)'a2{2) ) 
return 
end 
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•••*•  BLECTRONR 


•••••  electron  density  emd  gradient  for  Chapman  layer  model 


Input 

•  real's 

*  r  =  distance  from  Earth  centre 


*••**  output 

•  real's  • 

•  en  =  electron  density  (electrons  per  cubic  cm)  * 

•  dlffr  »  derivative  In  vertical  direction  • 


subroutines  re<iulred 
*  CHAP,  DCHAP  and  DET 


*****  cOTSDon  block  ***** 
*  DAT  (defined  In  HASEL) 


subroutine  electronr(r,en, dlffr) 

Implicit  real's  (a-h,o-z) 
dimension  al ( 3 ) , a2 ( 3 ) , a3 (3 ) ,b(3 ) 

common  /dat/foE,hiaE,yiaB,£oEl,  hmFl ,ymPl,  f oF2 ,hmF2,ymF2,  re, model 
If (r. It. re)  then 
eNoO.O 
dlffr-0.0 
else 

If (max(foE,foFl) .lt.l.d-7)  then 
xN»foP2**2/S0.6d-6 

eN»*N*chap( 1 .do, 1 . 4142d0* (r-ro-hmP2 ) /ymP2 ) 
dlf fr*l .4142*xN*dchap(l .dO, 1 . 4142d0* ( r-re-hmF2 ) /ymF2) /yinF2 
else 
al(l)-l. 

a2 ( 1 ) =ohap ( . 5d0 , 2 . dO* (hfflS-hmFl ) /ymPl ) 
a3 { 1 ) =chap { I . dO . 1 . 4 14 2dQ* (hmE-hmF2 ) /ymF2 ) 
al (2) =chap( . 5d0,2 .dO* (hmPl-hmE) /ymE) 
a2(2)»l. 

a3  ( 2 )  -chap  <  1  .dO ,  1 .4l42dQ*  (hmPl-hiaP2 )  /ymPa ) 
al  ( 3 ) -chap  ( .  SdO ,  2  .  dO*  (hBiK2-hmE) /ymE ) 
a2  ( 3 )  -chap( .  SdO ,  2  .dO*  (hjnF2-hmPl )  /ymPl ) 
a3(3)=l. 

b(l)  =  foE"2/80.6d-6 
b(2)=foPl"2/80.6d-6 
b(3)-foP2**2/80 .6d-S 
de-det ( al , a2 , a3 ) 
cl-det (b, a2 , a3 ) /de 
c2-det(al,b,a3) /de 
c3>det(al,a2,b) /de 

oN«cl*chap( .SdO, 2 .dO* (r-re-hmE) /ymE) 

&*c2*chap( .5dO,2.dO*(r-re-hmPl)/ymPl) 
t+c3*chap(l.d0, 1 .4142dQ* (r-re-hmP2 ) /ymP2 ) 

dlf fr«2 .*cl*dchap( . SdO, 2 -dO* (r-re-hmE) /ymE) /ymE 
t+2 .*c2*dchap( .SdO, 2 .dO* (r-re-hmPl) /ymFl ) /ymFl 
t+1.4142*c3*dchap(l.d0,1.4142d0'(r-re-hmP2)/ymP2)/ymP2 
endlf 
endlf 
return 
end 


""*  CHAP 


*****  chapman  layer  function 


*****  Input  •*•*••••••**••**••••••************** 

*  real*0 

*  C  -  1  In  P2  layer 

*  -  1/2’ In  B  and  PI  layers 

*  X  •  height  above  peak  In  units  of  scale  height 


*****  output  *•*•****•*•*•••••••**' 

•  real* 8 

*  chap  ■  value  of  Chapman  function 


real*8  function  chap(C,x) 
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implicit  real*8  (a-h, o-z) 
xx-exp( -X) 

ie(aba(xx) .gt.l.dSO)  than 
chap»0 . 0 
else 

ohap»exp(C» (l.dO-x-xx) ) 
endlf 
return 
end 

*****  DCHAP  ********•**********' 


*****  derivative  of  chapman  layer  ftmctlon 


*****  input  ******* **************************************** ************* 

*  real*8  • 

*  C  =  1  in  P2  layer  • 

*  *  1/2  in  E  2Uid  FI  layers  * 

*  X  •  height  above  peak  in  units  of  scale  height  * 


*****  output  •***•*••***************•••*' 

•  real*8 

*  dchap  *  derivative  of  Chapman  function 


real*e  function  dchap (C,x) 
implicit  real»8  (a-h,o-2) 
xx=exp( -X) 

if (ab8(xx) .gt.l.dlO)  then 
dchapsO . 0 
else 

dchap=C» (exp( -X) -l.dO ) *exp(C* (1 .dO-x-xx) ) 
endlf 
return 
end 

MAGNETIC  *•**********•***••*••***********« 


*****  magnetic  fields  and  derivatives  for  ao  offset  dipole  model 


*****  input  •*•*•**•*••*•**•**••••***•♦•••••***•••****•*••**••** 

*  real*8 

*  x{*)  =  Cartesian  coordinates  (x(3)  axis  in  the  Earths  axis  of 

*  rotation  and  x(l)  axis  in  the  plane  of  zero  longitude). 
*****  output  *************************************************** 

*  real*8 

*  v(*)  *  magnetic  field  components  in  Cartesian  system 

*  ryx(*,*)  =  derivatives  of  manetlc  field  comi>onents 

*  ha  *  magnitude  of  magnetic  field 

*  dip  *  magnetic  dip  angle 

*  del  «  declination  of  magnetic  field 


*****  common  block  •••*• 
•  MAG  (defined  in  HASEb) 


sub rout Ine  magnet ic { x , v , ryx , ha, dip, dc 1 ) 

implicit  real*8  (a-b,o-z) 

dimension  x(3 ) , v( 3) , ryx(3, 3 ) 

common  /mag/fhs, fr,xm,ym, zm,px,py.pz, Ir 

xmt-x(l ) -xm 

ymt«x(2)-ym 

zmt*x(3)-zm 

pm>>xmt*px*ymt*py+zmt*pz 
r  r  ■  xmt  •  xmt  ♦  ymt  •  ymt + zmt  *  zmt 
rasqrt(rr) 
rp«l./(rr*r) 

v(i)«(px-3.*^*xmt/rr)*rp 
V ( 2 ) - ( py- 3 . •pm*ymt / rr ) •  rp 
v{3)«(pz-3.»pm*zint/rr)*rp 
ra-sqrt(x(l)*x{l)+x(2)*x(2)*x(3)*x(3) ) 
ha»sqrt (v(l)»v(l) +v(2)*v(2) ♦v(3)»v(3) ) 
vr-v(l)*x(l)+v(2)*x(2)*v(3)*x<3) 

<lip>-45.*asln(vr/(ha*ra) )/atan(l. ) 

dcl-45.*atan( {v(2)*x(l) -v(l)«x{3) >/ (v(3)»ra-vr*x<3) ) ) /atan(l. ) 
ryx (I, 1 ) « (-3 . •tn/rt*6 . •pBi*xmt*x»t/ (rr*rr) -3 . •px*xmt/rr ) •rp 
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<t-3 .  *xsat*  ipx-3  .  •jan*xmt/rr )  ‘rp/rr 
ryx(l,  2 ) «  (-3 .  'py^xint/rr  +  e  .•iaii*ymt»xint/  (rr*rr)  )»rp 
s-3 . 'ymt*  <px-3 .•pm*xmt/rr) *rp/rr 
ryx (1, 3 )  =  (-3  .  *p2»xint/rr+6  . •poi»zmt*xiat/  (rr*rr )  ) *rp 
&-3 .  'znit*  {px-3  .  •pni*xmt/rr)  ‘rp/rr 
ryx  {2 , 2 ) «  ( -  3 .  ‘pm/rr+e  .  *pm*Yint*yint  /  ( rr*rr)  -  3  .  ‘py^ymt /rr)  *rp 
t-3.»ymt*(py-3  . •pin*ymt/rr) ‘rp/rr 
ryx(2, 3)»  (-3  .  *pz*yint/rr+6  .•pm»zint*ymt/  (rr*rr) )  ‘rp 
&-3.»zjnt*(py-3  . *jan*ymt/rr ) ‘rp/rr 
ryx{3, 3) »  (-3  .‘pm/rr+e  .*ian*zait»zmt/  (rr*rr)  -3 .  'pz^zmc/rr)  ‘rp 
S-3  .  “zint*  <pz-3  .  •txn*zmt/rr)  ‘rp/rr 
ryx(2,l)=ryx(l,2) 
ryx ( 3 , 2 ) «ryx ( 2 , 3 ) 
ryx(3.1)-ryxCl,3) 
return 
end 

PUNC 


****  right  hand  side  of  Haselgrove  equations 


Input  ••••••••••**• 

real*8 

voc(*)  »  value  of  solution  vactor  at  which  right  hand  required 
Integer 

n  =  number  of  components  In  vectors  vac  and  f 

output  •**•*••••*•*••****»**»»»*»•••***••*•***»*»••*••*»••< 

raal*0 

f(*)  •  right  hand  side  of  Haselgrove  equations 


*•**  subroutines  required 
ELECTRON,  MAGNETIC 


'•••  common  block  *•••• 
HAG  (defined  In  HASEL) 


*»*•  notes 

In  the  current  code  (HASEL),  y  has  nine  components.  In  order,  these 
are : 

the  Cartesian  coordinates  at  the  current  point, 
the  ccaoponents  of  the  wnve  vector, 
the  group  path, 
the  Doppler  shift, 
and  the  phase  path. 


subroutine  func (vec, f ,n) 

Implicit  real's  (a-h,o-z) 

dimension  f (n) ,vec(n) ,v(3) , ryx(3, 3) ,dnx(4 ) 
common  /mag/fhs,fr,xm,ym, zm,px,py,pz, Ir 
call  electron(vec,an,dnx) 
rxt»0 .06d-5*dnx(4 ) / (f r'fr) 
rx»8.06d-S*en/(fr»fr) 
vv«vec(4) '•2+vec(5) ••2+vec(6)**2 
If (ab8(fha) .gt.l.Oe-20)  then 
call  magnetic (vac, V, ryx, ha, dip, del) 
vu«  ( vec  ( 4 )  * V  ( 1 )  +vec  { 5 )  * v  ( 2 )  ♦vec  ( 6 )  •  v  ( 3 )  )  /  ha 
ry«ha 
ry2*ry*ry 
ga»-vu**2 
gga-ga/w 
al»l , -rx-ry2 
be» .5*ry* (1 .-ga) 
fac-sqrt(abs(be**2-al»ga) ) 

Ifdr.aq.l]  then 
rp--ga/ (be+fac) 
rj«4.*fac 
else 

rpa- (be+fac )/al 
rj»-4 .‘fac 
endif 
rp2»rp»*3 

rIt«-2  .*rx*vu*  (rp*ry-l . ) 
rl-ry*(l.-ry2)*rp2-2.«(l.-rx-ry2)*rp-ry 
rm»2.*rx*rp*(rp*ry-l. ) 
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do  l»l/3 

f  ( 1 )  »r  j *vec  ( i +3 )  - r)c»v  <  1 ) 
dxx>8 .06e-5*dnx(l) / (rr*fr) 
f  ( 1+3  )=rl»<ixx+  (rlt*v«c(4 )  +rm*v(l ) )  ‘ryxd,  1) 
fc+  (r)t*vec(5)  +nn*v(2)  )  ‘ryxC  2,  i)  ♦  (rlt*vee  (6)  ♦rm*v(3  ) )  •ryx(3,  i ) 
enddo 
rxl»l .-rx 
a»al+rx»ry2»gg 
b«rxl*al- . 5*rx*rY2* (1 . -gg) 
c-rxl*(rxl**2-ry2) 
amu»{b+ir*8grt (aba (b**2-a*c) ) ) /a 
amlxemu^a-b 
einsssgrt(abB(etBu) ) 

If (aba (ami) .gt .1 .Od-9)  then 
aa-1 . -1 . 5»rx-l . S»ry2+2 . •rx^rya'gg 
bb»rxl* (1.-3. ‘rx) -2 . •ry2  +  1 . 5*rx*ry2* ( 1 . +gg) 
cc--1.5*rx*rxl»*2-ry2*( .S-rx) 
emua  (aa  *eiBu  *ems  -bb*ama  +cc  /ams }  /eal 
ax«-l .+gg»ry2 
bx--2 .‘rxl* .5»ry2*(l. +gg) 
cx«-3 . *rxl**2+ry2 

emx» .25* ( -ax*em8**4+2 .•bx*em8**2-cx) / (a*ein8**3-b*amB) 
else 

emu-l ./sqrt (aba (1 . -rx) ) 
anxs- . 5*efflu 
endif 
else 

em8«8<irt(ab8(l.-rx) ) 
do  1>1,3 
f (i)-vac(i+3) 

dxx-8.06e-5*dnx(l)/(fr*fr) 
f  (1  +  3 )»- . 5*dxx 
enddo 

emual . /ema 
emx3-  .  5*emu 
endif 

coaa*  ( vec  (4 )  ‘f  ( 1 )  +vec  ( 5 )  •  f  ( 2 )  +vec  ( 6 )  *f  ( 3 ) )  / aqrt  ( w) 
f (7 ) semu'cosa 

f  (8  )  =  -co8a»fr*rxt»enix/3  .o5 
f (9 ) >ema*coaa 

if (abs(fh8) .gt.l.Oe-20)  than 
do  l=l<n 

f (l).real(lr)*f (1) 
enddo 
endif 
return 
end 


RKF 


Runge-Kutta-Feblberg  schema  to  solve  dy/dts  f 


Input  •****••••*•**•••••*••••••••••••••••••••***•*••• 

real*8 

t  s  Initial  value  of  solution  parameter  (changed  on  exit) 

y(*)  »  initial  value  of  solution  vector  (changed  on  exit) 

h  •  Initial  estimate  of  change  In  t 

tol  «  maxim  error  that  can  be  tolerated  In  this  step 

hmlns  minimum  h  that  can  be  tolerated 

hmax«  maximum  h  that  can  be  tolerated 

Integer 

n  a  number  of  components  In  vector  y 
output 
real*8 

t  M  exit  value  of  solution  parameter 

y(*)  a  estimate  of  solution  vector  at  exit  t 

b  a  estimate  of  b  for  next  step 


subroutine  required  ••*•*••*••*••••••••••••***•*••***••****•*•’ 

PONC 

This  calculates  the  value  of  f  corresponding  to  y.  There  are  three 
arguments.  The  first  Is  y  (real*8),  the  second  Is  f  (rsal*8)  and 
the  third  is  the  number  of  components  In  arrays  y  and  f  (Integer). 
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aubrout Ine  rkf ( t , y , n, h, tol , hmla, hmax) 

Implicit  raal*8  (a-h, o-z) 

dimension  al (11) , a2 (11) ,a3 (li) ,a4 (H) ,a5 (11) , a6 (11) ,y (n) ,yt (11 ) 
iendaO 

99  lends  iend-fl 

call  funo(y,al,n) 
alp=0 . 0 
do  Isl.n 
alp=alp+al ( i) ••2 
yt(l)=y(l)+.25*al (l)*h 
enddo 

alpsmax(sqrt (aba (alp) ) , 1 ■ Od-32) 
call  func(yt,a2,n) 
do  Isl.n 

yt(i)«y(l)+( .09375*al(l)».28125*a2(l) )*h 
enddo 

call  func (yt,a3,n) 
do  isl,n 

yt(l)=y(l)+( .87938097406d0»al(l)-3.2771961766d0*a2 (i) 
&+3.32089212563d0*a3(i) )*h 
enddo 

call  func (yt, a4,n) 
do  i«l,n 

yt(l)»y(l)+(2.0324074074dO«al(l)-8.»a2(l)+7  1734892787Sd0*a3 ( 1) 
S- .2058966a616d0«a4(i) )*h 
enddo 

call  func(yt,a5,n) 
do  isl/n 

yt(i)=y(l)+(-.29C29629629d0*al(i)+2.*a2(l)-1.38167641326d0«a3{l) 
<i+.4S297270955d0*a4(l)-.275»aS{l) )*h 
enddo 

call  func (yt, a6, n) 
erxsO .0 
do  i-1,3 

erz- . 277777777777778d-2»al(i) - . 299415204678363D-l»a3 ( i) 
fc- .291998936735779D-l»a4(l)+.2D-l*aS(l)+.36363o363636364D-l*a6(l) 
orxserx«-erz»*2 
enddo 

erxssqrt (erx) 
holdsh 

If (erx.gt.O.OOl'tol)  then 

as . 84»8qrt ( aqrt (aba (tol/erx) ) ) 

h«a*h 

else 

h=2 .‘h 

endlf 

if (h. It .hmln/alp)  hshmin/alp 
if (b.gt -hmax/alp)  hshmax/alp 
if (erx.gt.tol.and. lend.lt.8)  then 
go  to  99 
else 

do  lsl,n 

y(l)sy(l) t( .1185185ie52*al(l)+.51898635478d0*a3(l) 
s+.S0613149034dO*a4{l)-.18»a5(l)<-.03636363636*a6(l) )*hold 
enddo 
tst+hold 
return 
endlf 
end 
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•*•••  ELDEN 


•*•••  aeneral  electron  density 


*«••*  Input  **••***•*•**••••••*••••••*•••••••••*•••••***•**••••»• 

*  real»8 

*  t  »  time  from  start  (sec) 

*  elong  >  longitude  (deg)  of  sample  point 

*  elat  <  latitude  (deg)  of  sample  point 

*  r  3  height  above  Barth  centra  (Kms) 

* 

*  Note  that  hmE,hmFl  and  hmfS  should  be  set  to  suitable  values  to 

*  allow  the  labelling  of  rays 

**#»*  output  **************** ***•*••* ************ **************** 

*  real*e 

*  elden  «  electron  density  (electrons  per  cubic  cm) 


•»•••  common  block  (optional) 
*  DAT  (defined  In  HASEh} 


real*8  function  elden(t, elong, elat, r) 

Implicit  real*8  (a-h,o-x) 

camsion  /dat  /f  oE ,  hmE ,  ymE ,  f  oPl ,  hmPl ,  ymFl ,  f  oF2 ,  hniFS ,  ymF2 ,  re ,  mode  1 

foF2-12 . 

hmF2o350. 

yinP2»110. 

vela . 04 

x=(r-re-hiaP2-vel*t)  /yinF2 

eldea=(foP2**2/80.6d-6) *chap(l.d0,1.4142d0»x)+0.«olong»elat 

return 

end 


28 


UNCLASSIFIED 


UNCLASSIFIED 


SRL-OIJI-TR 


DISTRIBUTION 

Page  No. 

Defence  Science  and  Technology  Organisation 

Chief  Defence  Scientist  ) 

Central  Office  Executive  )  1 

Counsellor,  Defence  Science,  London  Cnt  Sht  Only 

Counseltar,  Defence  Science,  Washington  Cnt  Sht  Only 

Scientific  Adviser,  Defence  Central  1 

Navy  Office 

Navy  Scientific  Adviser  1 

Air  Office 

Air  Force  Scientific  Adviser  l 

Army  Office 

Scientific  Adviser,  Army  1 

Defence  Intelligence  Organisation 

Scientific  Adviser  1 

Surveillance  Research  Laboratory 

Director,  Surveillance  Research  Laboratory  1 

Chief  High  Frequency  Radar  Division  1 

Research  Leader,  Jindalee  Operational  Radar  Network  1 

Head,  Radar  Processing  and  Tracking  1 

Head,  Ionospheric  Effects  1 

Head,  Radar  Technology  and  Systems  1 

Head,  HF  Radar  Engineering  1 

Graf^ics  and  Oocumerrtation,  High  Frequency  Division  1 

Dr  J.A  Benn^,  Ionospheric  Effects  1 


UNCLASSIFIED 


29 


SRL-0131-TR 


UNCLASSIFIED 


DISTRIBUTION  (Contd) 

Chief  Microwave  Radar  Division  Qnt  Sht  Only 

Research  Leader  Microwave  Radar  1 

Research  Leader  Radar  Projects  ^ 

Author  ^ 

Electronics  Research  Laboratory 

Chief  Communications  Division  Cnj  Sht  Only 

Research  Leader  Military  Communications  1 

Research  Leader  Communications  Countermeasures  1 

Dr  K.J.W  Lynn,  Radiowave  Propagation  1 

Dr  R.H  Clarke.  Radiowave  Propagation  ^ 

Libraries  and  Information  Services 

Australian  Government  Publishing  Sen/ice  1 

Defence  Central  Library  -  Technical  Reports  Centre  1 

Manager,  Document  Exchange  Centre  (For  Retention)  1 

National  Technical  information  Service,  United  States  2 

Defence  Research  Information  Centre.  United  Kingdom  2 

Director  Scientific  Information  Services,  Canada  1 

Ministry  of  Defence,  New  Zealand  1 

National  Library  of  Australia  ^ 

Defence  Science  and  Technology  Organisation  Salisbury.  Research  Library  2 

British  Library  Document  Siqsply  Centre  ^ 

Library  Defence  Signals  Directorate,  Meboume  1 

Spares 

Defence  Science  and  Technology  Organisation  Salisbury.  Research  Library  6 


3« 


UNCLASSIFIED 


Page  Classtfication 
UNCLASSIFIED 


Department  of  Defence 


DOCUMENT  CONTROL  DATA  SHEET 


Privacy  MarkingfCaveat 
(of  Document) 
_ UIA _ 


1a.  AR  Number  1b.  Establishment  Number 
AR-008-171  SRL-0131-TR 


4.  Title 

A  GENERAL  PURPOSE  IONOSPHERIC 
RAY  TRACING  PROCEDURE 


2.  Document  Date 

August  1993 

3.  Task  Number 

DST  91/027 

5.  Security  Classification 

6.  No.  Of  Pages 

34 

m  nn  nn 

7.  No.  of  Refs. 

5 

Document  Title  Abstract 

S  (Secret)  C  (Confi )  R  (Rest)  U  (Unclass) 


8.  Author{s) 


'  For  UNCLASSIFIED  docs  with  a  secondary  distrtoution 
LIMITATION,  use  (L)  in  document  box. 


9.  Downgrading/Delimiting  Instructions 


C.J.  Coleman 


1 0a.  Corporate  Author  and  Address 

Surveillance  Research  Laboratory 
PO  Box  1500 
SALISBURY  SA  5108 


10b.  Task  Sponsor 


12.  Secondary  Distribution  of  this  Document 


N/A 


1 1 .  Officor/Position  responsible  for 


Security . N/A... 

Downgrading . N/A... 

Approval  for  Release . P.SRI^. 


APPROVED  FOR  PUBLIC  RELEASE 


Any  enquiries  outside  stated  limitations  should  be  referred  through  DSTIC,  Defence  Information  Services, 
Department  of  Defence,  Anzac  Park  West,  Canberra,  ACT  2600. 

13a.  Deliberate  Announcement 


No  limitation 

1 3b.  Casual  Announcement  (for  citation  in  other  documents) 

No  Limitation 

□ 

Ret.  by  Author ,  Doc  No.  and  date  only. 

14.  DEFTEST  Descriptors 

Ionospheric  propagation.  Ionospheric  properties,  Ray  tracing, 
Computer  ^plications 

1 5.  OISCAT  Subject  Codes 

0401,  1205 

16.  Abstract 

A  general  purpose  ionospheric  ray  uacing  procedure  is  developed.  The  procedure  is 
based  on  a  numerical  solution  to  the  Haselgrove  ray  tracing  equations  and  includes 
adrulations  of  ionospherically  induced  Doppler  shift.  The  procedure  can  handle  a  wide 
variety  of  ionospheric  descriptions  and  includes  the  effect  of  the  Earth's  magnetic  field. 


Page  Classification 

UNCLASSIFIED 


16.  Abstract  (CONT.) 


17.  Imprint 

Surveillance  Research  Laboratory 
PO  Box  1500 
SALISBURY  SA  5108 


18.  Document  Series  and  Number 

SRL-0131-TR 


19.  Cost  Code 


440/828309 


20.  Type  of  Report  and  Period  Covered 


TECHNICAL  REPORT 


El.  Computer  Programs  Used 


FORTRAN 


(22.  Establishment  Fite  Reference(s) 


N/A 


23.  Additional  information  (if  required) 


